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1. INTRODUCTION 

The continuum description of momentum and en- 
ergy transport in gases, based upon Newton- Stokes— 
Fourier constitutive relations, can become inaccurate in 
rarefied or highly nonequilbrium regimes, i.e., regimes 
in which the Knudsen number An (= A/L, where A is 
the gas mean free path and L is the characteristic system 
or gradient length) is no longer small. The Burnett 
equations, which represent the order— An 2 solution to 
the Boltzmann equation, ostensibly provide a means of 
extending continuum formulations into the transitional 
Knudsen regimes (~ An < 1). 

The accuracy and validity of the Burnett equations, 
however, have not been Firmly established. As has 
been noted by several authors, the asymptotic series 
expansion of the molecular distribution function — from 
which the Burnett equations are derived - has unknown 
convergence properties for finite An . 1,2 The Burnett 
equations can also lead to Second— law impossibilities, 
such as heat flux in an isothermal gas . 3 Furthermore, the 
Burnett equations increase the order of the differential 
equations that govern momentum and heat transport in 
the gas. Additional boundary condition information is 
required to fully close the problem - yet such informa- 
tion is generally not available from physical principles 
alone. 

Because of these issues, it is generally held that 
the Burnett equations arc valid only in regimes in which 
the Navier-Stokes-Fourier level of approximation al- 
ready provides an adequate description of transport, i.e., 
regimes in which the Burnett contributions represent 
a small perturbation to heat and momentum transport. 
Such conditions can be representative of high-Mach 
number flows, for which application of the Burnett 
equations appears to have been the most successful . 4 7 
On the other hand, there is not a broad understanding 
of the accuracy of the Burnett equations when applied 
to slow -moving, nonisothermal flow (SNIF) condi- 
tions. As noted by Kogan, ‘thermal stresses’ (fluid 
stresses resulting from temperature gradients - which 
are predicted by the Burnett equations) could become 
a significant convection mechanism in buoyancy -free, 
nonisothermal gases . 8,9 Indeed, it has been recently 
suggested that thermal stress convection could affect 
the growth of crystals in micro gravity physical vapor 
transport experiments. IOjI 1 

The work presented here consists of a theoreti- 


cal and numerical examination of thermally -induced 
stresses and flows in enclosed, highly nonisothermal 
gases under buoyancy— free conditions. A central objec- 
tive has been to identify a strategy in which stress and/or 
convection effects, as predicted by the Burnett equations, 
could be isolated and measured in microgravity -based 
experiments. Because of the questionable veracity of 
the Burnett equations, a second objective has been to 
test Burnett predictions of nonisothermal gas stress and 
convection with the exact description provided by the 
direct simulation Monte Carlo (DSMC) method. 

2. PREDICTION OF THERMAL STRESS 
CONVECTION 


The initial phase of the project was aimed at cal- 
culation, using continuum and DSMC methods, of gas 
convection in two dimensional nonuni formly heated 
rectangular enclosures. Typically, two adjacent surfaces 
of the enclosure were modeled as adiabatic, zero— stress 
surfaces (i.e., planes of symmetry), and the other two 
adjacent surfaces were maintained at specified temper- 
ature distributions with one surface transferring a net 
amount of heat to the gas, and the other transferring the 
heat from the gas. 

The continuum formulations of momentum and en- 
ergy transport arc identical to Navier-Stokes-Fourier 
models, with the exception of the Burnett stress tensor 
in the momentum equations and the creep and jump 
boundary conditions. For the conditions examined here 
(i.e., slow-moving flow, with Re^ 1 ) the only sig- 
nificant terms in the Burnett stress tensor relations will 
be those involving temperature gradients. This thermal 
stress component appears as 8,1 1 
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in which // is the dynamic viscosity, R is the gas 
constant, andu>* and 1 U 5 are dimensionless, order— unity 
coefficients which depend on the interaction potential of 
the molecules. The creep and jump boundary conditions 
appear 
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where n is the outward normal and the dimensionless co- 
efficients cs and ct depend on the thermal and momen- 
tum accommodation properties of the surface. 12 Nu- 
merical solution of the governing equations was accom- 
plished using the SIMPLER algorithm of Patankar. 13 
Coefficients corresponding to hard -sphere molecules, 
which gives a temperature- dependent viscosity of // ~ 
T 1 / 2 , were used in the computations. 

Direct simulation Monte Carlo calculations of 
hard— sphere gas convection and heat transfer were 
accomplished using the standard procedure developed 
by Bird. 14,15 The cell size was nominally set to 0.1A o , 
where A 0 represents the mean -free -path at the equi- 
librium state of the system, and 10 — 20 molecules were 
assigned per cell. Simulations were conducted for a 
Knudsen range of Kn — 0.01 — 0.2. Because thermal 
creep and stress flows will be on the order of Kn times 
the mean molecular velocity, resolution of the flows 
using DSMC required simulation times on the order of 
10 n -10 7 time steps. 

Our continuum and DSMC calculations to date in- 
dicate that it would be very difficult to create conditions 
in the enclosure that result in measurable thermal stress 
flows that are comparable to or larger than thermal 
creep flows, and simultaneously maintain the Kn <0.1 
regime required of the Burnett equations. With the ex- 
ception of the pure continuum limit (Kn -> 0, under 
which thermal stress vanishes), elimination of ther- 
mal creep cannot be accomplished by maintaining the 
hcated/cooled walls at uniform temperatures. Rather, 
the discontinuity (or jump) between the surface and ad- 
jacent gas temperatures — which will be proportional to 
Kn and the local normal temperature gradient - will 
lead to nonuniform gas temperatures along the nonuni- 
formly heated surfaces. For all realistic values of An, 
thermal creep flows generated by the temperature jump 
effects were substantially larger than those resulting 
from thermal stress. 

To minimize the effects of creep, we performed 
additional simulations in which the temperature dis- 
tributions along the heated/cooled surfaces were as- 
signed to provide, for a given An, nearly uniform gas 
temperature adjacent to the surface. Surface tempera- 
ture distributions were determined from solution of the 
gas conduction equation with uniform gas temperature 
boundary conditions along the heated/cooled surfaces, 
and subsequent application of the solution into Eq. (3) 
to predict T w (x). This approach imposed a surface 
temperature on the heated surface which increased to- 
wards the junction with the cooled surface, with an 
opposite trend along the cooled wall. The effect of 
this strategy resulted in thermal creep flows that were 
confined about the hot/cold junction, and left a bulk, 




Fig. 1 : DSMC (top) and continuum (bottom) results 

thermal— stress driven convection pattern in the bulk 
of the enclosure. However, this strategy would also 
lead to highly nonequilibrium conditions in the vicinity 
of the hot/cold surface junction, for which the Burnett 
equations are not expected to hold. 

The comparison between continuum/Bumetl and 
DSMC predictions of convective flows in the enclo- 
sure has been inconclusive. DSMC calculations have 
shown convective flows that are qualitatively similar 
than those predicted by the continuum model. As an 
example, we show in Fig. 1 plots of velocity vec- 
tors and isotherms, calculated using the DSMC (top) 
and continuum (bottom) models, for a Kn = 0.02 
hard- sphere gas contained in an enclosure with a rela- 
tive temperature difference on the hot and cold walls of 
2 (Tfj - Tc) /( Th + Tc) — 1. The heated wall is on the 
top, and the left and bottom walls are symmetry surfaces. 
A strip of the walls adjacent to the hot/cold junction, 
equal to 0.1 of the wall length, is held adiabatic, and the 
remainder of the walls have temperature distributions 
set to give isothermal gas conditions per the procedure 
discussed above. Velocity vectors corresponding to ther- 
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mal creep along the adiabatic surfaces (seen in the tight 
counterclockwise rotation in the upper right comer) have 
been removed from the plot to allow resolution of the 
thermal stress flows. The stress flow — as predicted by 
the continuum model — results in a clockwise rotation 
in the main body of the gas for the given conditions. 
A similar pattern is seen in the DSMC results. We 
cannot establish, however, whether the observed DSMC 
flows result from thermal stress, or are due to slip ef- 
fects at the walls. As mentioned above, the veracity of 
the Burnett relations is also questionable for the highly 
nonequilibrium conditions of the simulations. On the 
other hand, temperature profiles calculated via DSMC 
and continuum models show significant agreement. 

3. THERMAL STRESS IN I D HEAT 
TRANSFER 


A simpler situation in which to compare contin- 
uum/Bumett and DSMC predictions of thermal stress 
effects is offered by 1 — D heat transfer in a stationary 
gas. In this situation, the effects of thermal stress would 
be seen in the pressure distribution and normal stress in 
the gas. 16 

The computational domain was now taken to be 
a slab of gas contained between two parallel surfaces, 
separated by a distance A, with the surfaces at x = 0 
and L maintained at uniform temperatures of T C s and 
77/ s (with 77/s > Tcs\ respectively. In nondimen- 
sional form (with pressure and stress normalized with 
the equilibrium pressure P m and temperature by the 
equilibrium temperature T m ), the Burnett equation for 
the .r— directed, normal component of the stress tensor 
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In the above, 0 = PjP m% r* = r/P mi 0 = T/T m , the 
prime denotes differentiation with respect to £ = xj L, 
d — (4tt/ 3)(5/16) 2 = 0.4091, and the dimensionless 
u,* coefficients depend on the molecular interaction po- 
tential. Since the gas is stationary and buoyancy- free, 


the stress r will be a constant. In the limit of I\n — > 0, 
this gives the Navier- Stokes result of P = r = con- 
stant. For finite An, however, the additional source 
of thermal stress can act within the nonisothermal gas. 
The magnitude of the thermal stress will vary with po- 
sition - by virtue of the dependence of temperature and 
temperature gradient on position - and consequently 
pressure will vary to maintain a constant normal stress. 

The Burnett equations make no contribution to the 
heat flux for a stationary gas. Consequently, the gas 


temperature will be described by 

q* = constant = 0 l ^ 2 6 f (5) 


where q* is the dimensionless heat flux (= qL/kT m ). 
Equation (5) can be used to combine the last two terms 
in Eq. (4), which results in 
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where c 2 = (u; 3 - 2u? 5 )/2 - 0.9900 for hard-sphere 
molecules. 

Two separate effects - or regimes - on <j> can be 
anticipated from inspection of Eq. (C). One effect, which 
is discussed by Kogan 18 and Makashev 10 , derives from 
the fact that the derivatives of <f> will be multiplied by 
the small parameter (for near— continuum conditions) of 
Air. The solution to Eq. (G) could therefore exhibit 
‘boundary layers’ of width ~ Kn. It is shown below 
that this property, combined with appropriate boundary 
conditions, will allow for a limited description of the 
Knudsen layers adjacent to the surfaces. 

A second characteristic regime, as indicated by 
Eq. (G), would occur outside the Knudsen layers. By 
expanding <f> in a power series of Ah 2 , and neglecting 
all terms higher than Ah 2 , the pressure distribution in 
the bulk gas would be given approximately by 


<j> as r* + 


c\C 2 {Kn q’) 2 
T*0 


(7) 


As is evident from inspection of Eq. (7), thermal stress 
would create a pressure gradient in the gas, with pres- 
sure increasing towards the cooler regions in the gas. 
The gradient would be proportional to the square of 
Kn q* ~ y/0 dO/d{x/ A) - which can be interpreted as 
a Knudsen number based on the characteristic length of 
the temperature gradient (note that this quantity is inde- 
pendent of L). It should be emphasized that the effect 
predicted from Eq. (7) is fundamentally different than 
the pressure gradient created by ‘thermal transpiration’ 
of a gas in a tube with an imposed axial temperature 
gradient. 20,21 The latter is a result of thermal slip at the 
walls of the tube, and leads to a pressure that increases in 
the direction of increasing temperature. Thermal stress, 
on the other hand, results from the effect of temperature 
gradients on the molecular velocity distribution function 
within the gas. 

Although the thermal stress pressure gradient can 
be labeled ‘hydrostatic’ - since the gas is at rest - it 
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is distinctly different than that resulting from a gravita- 
tional acceleration in the .r- direction. Unlike the latter, 
thermal stress would not result in a difference between 
the normal forces acting on the hot and cold surfaces. In 
other words, the ‘pressure’ measured at the surfaces - 
which would physically represent the normal stress r — 
would be identical for both surfaces. 

Thermal stress, however, will result in a different 
value of the normal stress than that predicted from the 
Navier-Stokes level of approximation. This follows 
from conservation of energy requirements. In particular, 
the average pressure in the gas represents the equilib- 
rium pressure that would be attained if the walls were 
instantaneously made adiabatic. Since the equilibrium 
pressure is used to normalize the dimensional pressure 
P, this statement is equivalent to 

= 1 (8) 

Jo 

Regardless of the values of q* and Arc, the pressure 
distribution in the gas must satisfy the energy conser- 
vation constraint implied by Eq. (8). Consequently, the 
normal stress r* would be obtained as the eigenvalue 
to Eq. (6) such that the solution (for specified boundary 
conditions) satisfies Eq. (8). In general, this value will 
be different than the Navier-Stokes result of r* — 1. 

An approximate value for r* can be obtained by 
neglecting the effects of the Knudsen layers at the 
surfaces, for which the pressure distribution would be 
given by Eq. (7). To order An 2 , this gives 16 

r* « 1 - c,c 2 (Knq*) 2 (9) 

This relatively -simple approximation indicates that 
thermal stress will lower the normal stress in a closed 
system relative to that predicted from the Navier-Stokes 
level - although we note again that the effects of Knud- 
sen layers have been neglected in the analysis. 

The final elements required to close the prob- 
lem arc the boundary conditions for pressure. As is 
the case with the Navier-Stokes approximation, the 
boundary conditions for the Burnett equations should 
represent an extrapolation of the solution across the 
region, adjacent to the wall, where the solution is no 
longer valid. Makashev 19 and Schamberg 22 have pro- 
posed boundary conditions that are consistent with the 
order- K n 2 accuracy of the Burnett equations. The 
accuracy of these approaches, however, has not been 
well established. 23 Alternatively, order- An relations 
can be derived for the pressure ‘slip’ adjacent to a heated 
or cooled surface. 12,21 However, our work at this stage 
is primarily concerned with determining whether there 
are boundary conditions which, when coupled to the 



Figs. 2-4: DSMC and continuum pressure distributions 

Burnett equations, can reproduce DSMC predictions of 
pressure distributions in a nonisothermal gas. Therefore, 
the pressures at the hot and cold surfaces were taken 
to be parameters, and were chosen to provide the best 
agreement between theory and DSMC results. The ob- 
vious choice for the pressure at the surfaces will be the 
values determined from DSMC predictions. 

Comparisons of Burnett (via numerical solution of 
Eq. (6)) and DSMC predictions of pressure distribution 
appear in Figs. 2—4. Each plot shows P/ r = <P/r* vs. 
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x i — xjL for a fixed value of q *, with Ah a parame- 
ter. Two theoretical predictions of P/r are shown for 
each set of DSMC results. The first corresponds to the 
numerical solution of Eq. (6), with boundary conditions 
obtained from extrapolation of the DSMC -derived pres- 
sures to the surfaces and temperatures predicted from 
Eq. (5). The second represents the bulk gas thermal stress 
pressure distribution predicted from Eq. (7), which does 
not account for boundary effects. This latter predic- 
tion has been shifted by a constant to match with the 
full-Bumett solution at £ = 0.5. 

As is evident from the results, the pressure profiles 
show distinct Knudscn layers at both the cooled and 
heated surfaces. The drop in pressure at the cooled 
surface, and the increase in pressure at the hot surface, 
are both consistent with the predictions of pressure slip 
relations. 21 The pressure drop at the cold surface can 
be considerable for the conditions examined here — 
amounting to around an 8% decrease for An = 0.2 and 
q* = 1.5. 

The solutions of the Burnett equation, with DSMC 
derived boundary conditions, are seen to capture the 
essential features of the DSMC pressure distribution. 
In particular, the solutions provide a good description 
of the width and form of the Knudscn layers and the 
pressure distribution outside the layers. The difference 
between the theoretical and DSMC results is greatest at 
the edge of the cold -surface Knudscn layers, for which 
the theoretical model tends to overprcdict the pressure. 
This is most evident for the results corresponding to 
q* — 2.0 in Fig. 4. Nevertheless, the fact that the 
Burnett equations can resolve, to a reasonable accuracy, 
the Knudscn layers at the surface is somewhat surprising 
- especially when considering that the theory is based on 
the order— An continuum temperature profile. We also 
examined solutions to Eq. (G) using boundary values of (f> 
that were different than the DSMC results, and found that 
the exact, DSMC-derived boundary conditions provide 
the best overall agreement between Burnett equation 
predictions and DSMC results. 

The DSMC results for q* =2.0 appear to show 
a pressure distribution in the bulk gas that is described 
by Eq. (7). On the other hand, the pressure distribution 
for q* = 1.5 and Ah = 0.2 (Fig. 3) - for which 
Eq. (7) predicts a greater effect — is dominated by 
the Knudscn layers extending from the surfaces. To 
eliminate the effects of the Knudsen layers at the hot 
wall, we performed additional DSMC calculations in 
which the velocities of the incoming molecules at the 
hot boundary were sampled from the Chapman -Enskog 
distribution function for the fixed values of q* and An. 
By doing so, the hot surface now approximated an open 
boundary. The corresponding DSMC results showed a 



Fig. 5: DSMC and continuum normal stress 

pressure distribution in the bulk gas that was accurately 
represented by Eq. (7). 

A final comparison of theory and experiment can 
be obtained from the dimensionless normal stress. The 
simplified model of Eq. (9) indicates that t* should be 
a function primarily of Kn q*. Accordingly, we plot in 
Fig. 5 the DSMC values of r* vs. Kn q* for the seven 
different combinations of Kn and q* that were used 
in the closed system calculations (results of Figs. 2-4). 
Theoretical results correspond to the derived eigenvalues 
of Eq. (G) for the DSMC-derived boundary conditions, 
and to the approximation given by Eq. (9). 

The first point to make is that the predictions of r* 
from full solution of the Burnett differential equation are 
nearly equivalent to those obtained from the bulk -gas 
approximation of Eq. (9). Evidently, the decrease in 
pressure at the cold surface is compensated by the 
increase at the hot, so that the Knudsen layers have 
a small effect on the averaged pressure in the gas. 
Secondly, the primary dependence of r* on Ah q* is 
supported in the DSMC results at Kn q* = 0.1 and 0.2 
— which each correspond to two combinations of An 
and q*. As observed, the results are nearly identical 
at these points. Finally, the theoretical predictions 
are in excellent agreement with the DSMC results for 
Ah q* < 0.15, beyond which the theory overpredicts 
the decrease in r*. As can be seen from the results, the 
relative decrease in normal stress on the surfaces is quite 
small, i.e., r* = 0.975 for q* — 2.0 and An = 0.2, or 
a 2.5% decrease in ‘measured’ pressure at the surface. 
We should emphasize, however, that this decrease is 
still significantly larger than the numerical precision of 
the DSMC simulations. 

SUMMARY 

The project has sought to ascertain the veracity of 
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the Burnett relations, as applied to slow moving, highly 
nonisothermal gases, by comparison of convection and 
stress predictions with those generated by the DSMC 
method. The Burnett equations were found to provide 
reasonable descriptions of the pressure distribution and 
normal stress in stationary gases with a 1 -D temperature 
gradient. Continuum/Bumett predictions of thermal 
stress convection in 2— D heated enclosures, however, 
are not quantitatively supported by DSMC results. For 
such situations, it appears that thermal creep flows, 
generated at the boundaries of the enclosure, will be 
significantly larger than the flows resulting from thermal 
stress in the gas. 
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